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Scope of. Extended Contract NAS2-4151 


Work under Contract NAS2-4151 started on February 1, 1967, 
Results obtained through August, 1967 are summarized in 
"Phase I Report under Contract NAS2-4151" of September, 1967. 
Subject contract was extended through July, 1968. The main 
purpose of the extended contract was to check the validity of 
the approximate digital method for computing the response of 
blade flapping to random inputs, tentatively suggested in 
Phase X Report, by' comparison with NASA conducted simulator 
studies, to develop alternate methods if required and to 
extend the analysis to higher rotor advance ratios. This 
report summarizes the results obtained since September, 1967 
through July, 1968, during which period 12,9 man-months were 
expended . 



-i- 


Goncepts for a Theoretical and Experimental 
Study of Lifting Rotor Random Loads 
and Vibrations 
(Phase II) 


by Kurt H. Hohenemser 

and Gopal H. Gaonkar 

Washington University, St. Louis, Missouri 


Abstract : 

A comparison with NASA conducted simulator studies has shown 
that the approximate digital method for computing rotor blade 
flapping responses to random Inputs, tentatively suggested in 
Phase I Report, gives with Increasing rotor advance ratio 
the wrong trend. Consequently, three alternative methods of 
solution have been considered and are described in this report. 
An approximate method based on the functional relation between 
input and output double frequency spectra, a numerical method 
based on the system responses to deterministic Inputs and a 
perturbation approach. Among these the perturbation method 
requires the least amount of computation- and has been 
developed in two forms - the first form to obtain the response 
correlation function and the second for the time averaged 
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spectra of flapping oscillations. The range of validity of 
the first form has been ascertained by a comparison between 
the Runga-Kutta and perturbation response values to harmonic 
inputs and that of the second form by comparing the time averaged 
response spectra values obtained from the perturbation method 
and the NASA conducted simulator results. Such comparisons 
indicate that the perturbation scheme should provide reasonable 
approximations up to a rotor advance ratio of one at a Lock 
blade inertia number of four. 
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Notation 
|l x = E [x] 


E [f(x)] 


x, = x 1 (t 1 ),x 2 = x 2 (tg) 


Expected value of sample x 

Expected value of sample function f(x) 

Values of sample function x(t) at 
times t^ and t 2 respectively 

Time 


t = t 2 - t. 


t = 


t 1 + 1 2 


Time difference 


Average time 


Rxy^J.tg) 


f U 

f “ TW 
Af 


Cross-correlation function between 
sample time functions x(t^) and y(t 2 ) 

Frequency 

Frequency Interval 


X(f) 


S xy^ f 1' f 2^ 


Fourier transform of sample 
function x(t) 

E [x*(r 1 )Y(f 2 )] Cross-correlation function 

between sample frequency functions 
x*(f 1 ) and y(f 2 ), also called 
power spectral density 


h( T) 


Unit Impulse response function 


H(f) 


Frequency response function 
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P 

So 

H 

0 

a 

V f) 

(f) 

A(t) 
y( a> >t) 
y A ( w ,t) 

y c ( u) >t) 
y g ( u,t) 

y Ao (w,t) 

y ka {u,t) 


Modulating frequency 
Rotor angular velocity 
Advance ratio 

Blade flapping angle, positive up 

Mean blade angle of attack 

Power spectral density of mean 
angle of attack 

Average power spectral density of 
blade flapping angle 

Right hand side deterministic function 

i fil t 

Response of the system to the input e 

Response of the system to the 
input A(t)e ia/t 

Real part of y( cj ,t) 

Imaginary part of y(6>,t) 

Real part of y A (<J,t) 


Imaginary part of yCt^t) 
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c 1' c 2 ,d 1 and d 2 


a Q ,a^ j • • • > ^■j ^^2 * « • 

B(f) 

A(f) 

Superscripts: 

* 


Constants associated with the 
left hand side of the blade 
flapping equation 

Constants associated with A(t) 

Fourier transform of pit) 

Fourier transform of a (t) 


Conjugate complex 
Time differentiation 


Time average 
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1 . Introduction 

Of the various complications encountered when trying 
to apply to lifting rotors the stochastic methods developed 
to analyze airplane responses to atmospheric turbulence, 
we are concerned in this Phase II Report only with the time 
varying character of the system and of the random input. 

Though the theory Is developed in a more general form 
applicable to the response of time varying linear systems to 
certain types of non- stationary random inputs, the application 
is to the flapping response of rigid blades hinged to a rigidly 
supported hub. The blades represent in forward flight an 
approximately linear system with time variable periodic stiffness 
and damping. Because of the periodically varying relative 
flow velocity occurring in forward flight of the lifting rotor, 
the aerodynamic excitation of the blades cannot be represented 
by a stationary random process as in the case of frozen wing 
aircraft, but must be described as a non- stationary stochastic 

input . 

Non- stationary random inputs have been analyzed for a 
few engineering applications, for example for the response of 
airframes to random runway disturbances during decelerations 
after touchdown, Ref. (1), for the description of strong motion 
random earthquake excitation. Ref. (2), and for the response of 
spacecraft to time varying random excitations during the 
launching phase, Ref. (3). In these applications the system 
had constant parameters and could be represented by a time 
invariable transfer function, A flapping blade of a lifting 
rotor, however, because of the time variable periodic parameters, 
cannot be represented by such a time invariable transfer function. 

The general theory of non- stationary stochastic processes 
has been well established' as a direct extension of the corres- 
ponding theory of stationary random processes. References (4), 

(5), (6), (7). Rigorous solutions of responses to non- stationary 
random inputs thus far available are, however, restricted to 
constant parameter systems. References (6), (8). The complexity 
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of the analysis is due to the fact that, except in very special 
cases, closed form solutions valid over large time intervals 
do not as yet exist for differential equations with variable 
coefficients. When it is possible to find a rigorous solution, 
as in the case of the Bessel differential equation, the 
quadrature operations required to obtain for example the 
response auto-correlation function are quite involved even for 
stationary random inputs, Ref. (8). 

Since our literature survey has not uncovered prior 
work toward solving the response of a linear system with time 
variable parameters under non-stationary stochastic loading, 
an approximate method for moderate advance ratio was suggested 
in Phase I Report according to which it was assumed that 
both the excitation and the response can be considered to be a 
stationary random process modulated by a deterministic time 
function. It was further assumed that the general equation 
between the two- frequency input and output power spectral 
densities could be approximately solved as far as the time 
averaged single frequency response spectrum is concerned, 
by ignoring the relations between the diagonal and off-diagonal 
terms of the two- frequency power spectra and by considering 
only the relation between the diagonal terms of the input and 
response spectra. 

Since submitting Phase I Report consistent data have 
been received from the Simulator Computer Systems Branch 
of the NASA Ames Research Center, where the problem of random 
blade flapping response at various advance ratios had been 
simulated upon our request. When comparing the time averaged 
single frequency power spectra from the simulator study with 
the equivalent data from the approximate digital method suggested 
in Phase I Report, it was 'found that this method resulted with 
increasing rotor advance ratio in the wrong trend, so that the 
method cannot be accepted as an approximation. © 
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Ef forts were then directed toward solving the general 
equation between the two- frequency input and output power 
spectral densities including the relations between diagonal 
and off-diagonal terms. However, even for rather crude dis- 
cretization of this functional equation one must obtain a 
computer solution for many hundred simultaneous linear equations 
between complex variables. It was found that a first computer 
program established for this purpose did not yield convergence 
of the iterations and this attempt was then suspended, though 
further work on the computer program might still lead to a 
success. 

Next it was considered to obtain a solution based 
on the response to deterministic inputs .assuming zero initial 
displacement and rate of displacement. Once such response 
time histories have been obtained for a sufficient number of 
frequency intervals, it is a simple matter of a frequency 
quadrature to obtain time variable response mean square values. 

It was finally decided, before attempting an entirely 
numerical solution, to develop a perturbation method of solution 
which is an approximate analytical method for cases where 
the time varying parameters in the differential equations do 
not differ very much from their time mean. A numerical solution 
of the deterministic response problem with its inherent 
computer costs. can then be avoided. 

Much of the reliability of the perturbation method 
within its range of applicability stems from the fact that it 
repeatedly deals with constant parameter systems for different 
known input s . 

Numerical evaluations in this report are concerned 
with the determination of the range of validity of the 
perturbation method and its application to the problem of 
random rigid blade flapping. 
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2 . General Relations and Concepts for Non- Stationary Random 
Processes 

Non- stationary random processes can be represented 
either by a double frequency power spectrum or in some cases 
by an instantaneous time varying spectrum. In the first two 
sub -sections these two important concepts are briefly discussed. 

The third sub-section deals with a particular non- stationary 
random process obtained by modulation of a stationary process 
with a deterministic time function. In the subsequent section 
on random response analyses it will be assumed that the non- 
stationary input is of this form. 

2.1 Double Frequency Spectra 

The general relations between correlation functions 
and power spectral densities have been discussed in Phase I 
Report Section 2.1. For two non- stationary random processes 
with sample functions x(t), y(t) having zero means and having 
the sample Fourier transforms X(f), Y(f) the cross-correlation 
function is given by 

«xy< t V t 2 ) = E [ x ( t 1 )y(tg) ] = fj S X y( f V f 2> e ' 12lr(fltl " 

— 00 

2.1 

and the cross-power spectral density is given by the Inverse of 
eqn. 2.1 

S xy (f r f 2 ) = E [x*(f 1 )y(f 2 )J = JJ R xy (t 1 ,t 2 )e 12,r(f 1 fc 1 - 

—00 

2.2 

For a single non-stationary random process with sample function x(t), 
autocorrelation function ^(t-ptg) and power spectral density 
S x (f^,f 2 ) are related by 

R x (t r t 2 ) = E = Jj S x (f 1 ,f 2 )e -i2lr(f 1 t 1 " 

—do 


2.3 
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and its inverse 

ct> 

S x (f 1 ,f 2 ) = E [x*(f,)X(f 2 j] = jf R x (t 1 ,t 2 )e i 2 lr(f 1 t 1 ‘ 

2.4 


The spectral function S (f^fg) is in general for any combination 
of frequencies 2 a complex number and not physically realizable. 
For weakly stationary random processes 


VVtg) « R^tg - t,) 


2 .5 

2.6 


where $ (•••••) is the Dirac delta function with properties 


<$(t) « 0 

i(t) * oo 
00 

J6(t)at 


if t # o 
if t = 0 


2.7 


— 03 


j d(t)dt = 


1 , €>0 


so that for a function 0(t): 


A* - t 0 ) *(t)dt = *(t 0 ) 2.8 


Inserting 2.5 to 2.7 in eqn. 2.3 and 2.4 one obtains with 
t 2 - t-j =r and f| = f 2 * f 


R x ( T ) = E [klt^xttg)] 
S x (f) = E [x*(f)x(f)] 


= f 5 x ( f )e i2irfT df 

Jn x ( t ) e _i2lrfT d r 


2.9 


2.10 


Since X*(f)X(f) is real, the spectral function S (f) is also 

A 

real and physically realizable. 
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2.2 Instantaneous and Time Averaged Spectra 

In Phase I Report a specific example of a time averaged 
power spectrum was given in Section 2.3. Here the general 
concepts of instantaneous and time averaged power spectra will 
be discussed. In the definition of the double frequency power 
spectral density, eqn. 2.4, it is assumed that the 'sample 
functions x(t) are defined over an infinite time interval, 
though in actuality one has available only sample records 
observed over a finite duration. Furthermore, the double 
frequency spectrum is not a physically realizable quantity. 

The concept of an instantaneous power spectrum. Ref. ( 9 ), 
defined over a finite time interval and then time averaged, 
resulting in a physically realizable single frequency spectrum, 
allows to establish a relation between field observations and 
the theory of non- stationary random processes. 

The Instantaneous power spectrum S x (f,t) is, according 
to Ref. (9 ), defined as 

S x (f,t) = J R x (t - § , t + §)e " 12lrfT dr 2.11 

— ao 

It is also not physically realizable but it leads to a practical 
way of treating non- stationary random data by averaging over a 
sufficiently long time interval T . Denoting with R x ( r ) and 
S’ (f) the time averaged autocorrelation function and power- 
spectral density respectively, one obtains from eqn. 2 . 11 : 

S x (f) = ^ 5 T J s x (f,t)dt - J R x (T)e _12irfT clT 

2.1a 

where 

R x ( T > = £?« t / ^(t - | , t + f)dt 

— <“ 

- f S x (f)e 12,rfT 

-03 


df 


2.13 
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In order to relate the double frequency spectrum S^f-pfg) to the 
time averaged spectrum S (f) we substitute in eqn. °,3 t, «* t - 3- , 
t^ = t + ^ , so that 




oc 12Trt(f 2 - f,) iirr(f 2 + f,) 

(t - | , t + p = 7_/s x ( fl ,f 2 )e 


df 1 df 2 


—co 


Therefore from eqn. 2.13 

T/2 


V T ) = 


lim 

T-><® 


i2rt(f 2 - f,) iirr (f 2 + f,) 


df 1 df 2 dt 


lim 

T-^co 


lirr(f 2 + f 1 ) i y z X2irt(f 2 - f , ) 


f/s x (f,,r 2 )e " ' 


Urn (I rp p . lirT(f 2 + f 1 } sin(f 2 - f,)^ 

~ _// s x( f r f 2 ) e (f 2 - df i 


dtdf 1 df 2 


df 


2 

2.14 


In accordance with Ref. (7), p. 450 it is now assumed 
that the double frequency spectrum of the non stationary random 
process can be expressed as a sum of regular and singular masses: 


s x (f r f 2 ) ■ MVV + S s ( f i> a ( f 2 ' f i } 

where S (f^^) has no line masses on the line f^ = fg. Inserting 
this expression into eqn. 2.14 and noting that 


sin(f - fJirT 1 for f, = f 
T— ■» (f 2 - f , 1 )7* 0 for f 1 / f £ 


f, + f 2 

one obtains with: n = f 


M t ) = 


00 

■J- 


12irr f 


Sjf)e 


df 


and by comparison with eqn. 2.13 


S x (f) = S s (f) 


2.15 



- 8 - 


This important theorem says that the time averaged autocorrelation 
function R^( T ) and the time averaged power- spectral density 
S (f) are uniquely defined by the line masses of the double 

A 

frequency spectrum along the diagonal f-j » f^. If these line 
masses are zero, then S^f) = R^r) “ 0* In the following 
section it is shown that such line masses occur when the non- 
stationary process can be represented by a stationary random 
process modulated by a periodic time function. 

2.3 Stationary Random Processes Modulated by a Periodic Time 
Function 

Consider a sample function from a non- stationary random 

process 

z(t) = A(t)x(t) 2.16 


where A(t) is a deterministic periodic time function and x(t) 
a sample function from a stationary random process. 

V/e write the Fourier series for A(t) with the basic frequency f 
in the form 


A(t) 



ik2irf t 
o 


2.17 


A(t) must be real since it represents a physical quantity, so that 
c k “ c -k and we ' can write e 3 n * 2 * 17 also in the form 


A(t) 



-ik2irf t 
o 


Applying the first part of eqn. 2.3, one obtains for the auto- 
correlation function 

i2'rrf 0 (-kt 1 + ltg) 

2.16 

Introducing as before t^ « t - r /2, t 2 * t + r/2 and considering 

that V*' T ) is independent of t, since x.(t) Is the sample function 


R z (ti,t 2 ) - 
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of a stationary random process, eqn. 2.18 assumes the form 



* ±2irf 0 (1 - k)t + (1 + k)g 
°k c 1 e 

2.19 


Time averaging the autocorrelation function according to eqn. 2.13 
and considering 

lim . f i2Tf„(l - k)t 1 for 1 = k 


CD 


1 f i ^7rf 0 (l 

*-f/2 


dt = 


0 for 1 * k 


one obtains 


oo * i2*rrf kr 

R,(r) = R X (T) = E, c k c k e 


2.20 


From the first part of eqn. 2.4 one obtains for the power spectral 
density 


C J C P S x \ 


(h.1 


f - f ( 1 
2 o VJ 


+ p)\ 


a(f 2 - f, - f G 


2.21 


The time averaged spectral density, according to the theorem 
proven in section 2.2, is equal to the line mass of the double 
frequency spectrum along the diagonal = f^* so that with J = P 


S z (f) c > j S x (f - Jf e ) 2.22 

The same result can also be derived by inserting eqn. 2.20 into 
eqn. 2.12. 


X. 

(-j+p) 


and 
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3 . Response of Time Varying Linear Systems 

In Section 2.4 of Phase I Report the response auto- 
correlation function and the response power spectral density 
were given In terms of the time variable impulse response function 
and of the time variable frequency response function. It was noted 
that numerical solutions would be very difficult since in general 
neither the impulse response function nor the frequency response 
function are given analytically. 

Following some remarks by Sveshnikov in Ref. (8), p. 135 
the problem is here reformulated by introducing the response y(f,t) 
of the linear system to an excitation e ^ where the system is 
assumed to be at rest in Its equilibrium position at the origin 
of time. It Is shown how input-output relations between correlation 
and spectral functions can be expressed In terms of the particular 
solution y( f , t) . 


3.1 General Non- Stationary Random Input 

The response y(t) of a time varying linear system with 
an infinite operating time is given by 


y(t) 



h( T *t) X 


( t - r ) dr 


3.1 


Physical realizability requires that the impulse response function 


h( r ,t). = 0 for T <0. 


Stability of the system requires that 


CO 

fY 

— oo 


h( T ,t) 


d T < oo . 


The input x(t) applied to the system at time t -t can be either 
a deterministic function or a sample function of a stochastic 
process. 

Assuming a harmonic input 


i2irft 

x(t) = e 


3.2 
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the output, according to eqn. 3.1 is 

CO 


r i2irf(t - T ) 

y(f,t) = J h( T ,t)e dT 


-CD 


which can also be written as 

12irft 

y(f,t) = H{f,t)e 


with 


/-® -i27rfr 

H(f,t) = J h( T ,t)e dT 

CO 


3.3 


3.4 

3.5 


For a general input x(t) we write the Fourier transform inverse 

i2irft 

e X(f)df 3.6 

Since the response to the input e^^^ is y(f,t), see equation 
3.2 and 3.3, one can write the response to the input x(t) in the 
form 

y(t) = / y(f ,t)X(f )df 3.7 

-CO 

x(t) and y(t) are physical qualities and, therefore, real. From 
eqn. 3.6 it then follows that X(f) » X*(-f) and from eqn. 3.7 that 
y(f,t) = y*(-f,t). The latter equation can now also be written as 

y(t) - f y*(f,t)x*(f)df 
-00 

Defining correlation functions and power spectral densities 
according to eqns. 2.1 to 2.4 and inserting eqns. 3.6 and 3.7 
one obtains: 

00 

R (ti,t 2 ) = _//'y*(f 3 ,t 1 )y(f 4 ,t 2 )s x (f 3 ,f 4 )df 3 df 4 3.8 

^ -CD 

oo -127rf t 

R xy (t 1 ,t 2 ) = ff y(f 4 ,t 2 ) e 3 1 S x (f 3 ,f 4 )df 3 df 4 3.9 
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r°r S r° 12irf i t i , 

//S x (f 3 ,f 4 ) 1 J e V 

— CO —CO 

{ /° - 127 rf 0 t 

Lo/ e y( f 4»t 2 )at 


-i2irf p t f 

e 2 y(f 4 ,t 2 )dt 2 ^ df 3 df 4 


3.10 


r ( r° -i2irf P t ( 

^ ^(f-j ,f^) | J s ^ df^ 


3.11 


The time dependent mean square response, which gives in many 
applications adequate information, is obtained by setting in 
eqn. 3.8 t-j = tg = t. Thus 

OO 

Oy (t) = Ry(t,t) = J'J* y ( f )y( ^ 3 * ^ 4 ) ^^ 3^4 

— OO 

3.12 


If the Input Is a stationary random process the double integrals 
of the frequencies reduce to single integrals and one obtains 
for example for the response autocorrelation function 
00 

R y (t 1 ,t 2 ) = J y*(f,t 1 )y(f,t 2 )S x (f)df 3.8a 

— OO 

and for the time variable mean square 

OO 

®y 2(t) = V t,fc) = f y*( f ,t) y (f,t)s x (f)df 3 . 12 a 

—CO 

For constant parameter systems the frequency response function H^f) is 
independent of time, so that 

127Tft 

y(f,t) = H(f)e 


3.13 
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The mean square response is then given by 


<r y 2 (t) = Ry(t,t) 



- 12¥t(f, 

f 3 )H(f 4 )e 



S x ( ^3^ f 4 )df 3 df^ 

3.14 


and is time variable. For stationary random input this reduces 
to the time invariable mean square response 

CO 

<^ 2 = Ry(0) = yH*(f)H(f)S x (f)df 3.15 

— oo 

For time variable systems the particular solution 
y(f,t) has to be usually evaluated numerically. It is, therefore, 
preferable to avoid complex arithmetic in the numerical algorithm. 
Since the system is linear, the response can be expressed as 


y(f,t) = y 0 (f,t) + iy s (f,t) 3.16 

where y ff,t) and y 0 (f,t) are the responses to the inputs cos2Trft 
c s 

and sin27rft respectively. By expressing the Fourier transform 
of a sample function by 

X(f) = X 1 (f) + iX 2 (f) 

and then applying the definition 2.4 

one can easily show that the power spectral density can be expressed 
by 

■yv f 2 > = s xR ( f i’ f 2 ) +is xi< f v f 2 > 3 - 17 


with properties 




(f v f 2 ) = s xR (-f r -f 2 ) 

{ f 1 ,f 2 ^ = _S xI ( ~ f 1 ,_f 2 


) 


3.18 


3.19 
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VW - ff {y 0 ( f r t i )y o (f 2' t 2 ) 

“00 

+ y s (fvti)y B (f 2> t 2 )( s xR (f v f 2 )d fl df 2 
00 

ff j 

CO 

" y c^ f 1 ,t: 1 ^ y s^2 jt 2^ } S xi( f V f 2^ df 1 Cif 2 


+ 

— Co 


RxyCt^tg) = ff { y 0 (f 2 ,t 2 )cos 8rf 1 t 1 

— 00 

+ y s^ f 2 ,t 2^ Sin 27rf i } S xR^ f 1 ,f 2^ df 1 df 2 
.00 

ff { y 0 (f 2 »t 2 )sin 2 1 rf 1 t 1 
00 

- y s (f 2 ,t 2 )cos 2irf -j } S xl^ f 1 ' f 2^ df 1 df i 


+ 

» 

—00 


For stationary random input equation 3.20 reduces 
integral 

r° 

R y ( t 1 .t 2 ) = J { y c (f,t 1 )y c (f,t 2 ) + y s (f,t 1 )y s (f 


The time variable mean square is then given by 
cr y 2 (t) = R y (t,t) = J |y c 2 (f,t) + y s 2 (f,t)| s. 


— 00 


3.20 


3.21 


to the single 


,t 2 )}s x ( f ) d f 

3.22 


(f)df 

3.23 


an expression derived by Sveshnikov in Ref. (8), p. 136. 
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3.2 Modulated Stationary Random Input 

It is now assumed that the input can be represented 
by eqn. 2.16 and 2.17 so that the input is a stationary random 
process modulated by a periodic time function. In section 2.3 
it had been shown that such a random process leads to a special 
kind of double frequency power spectrum with line masses along 
the diagonal f-j => fg, see eqn. 2.21, so that an average power 
spectral density and an average autocorrelation function different 
from zero exist. Substituting the input spectrum eqn. 2.21 
into eqn. 3,8 for the response autocorrelation function, one 
obtains 



Jj v ( f i )y( f 2 ,t 2^ k 

— OO -K’ 


CD OQ 

£ £ 

-°°l= -» 


C, C 

k 


,4 


1 + f 2 


-f 0 (k + 1) 
1 



«(f 2 - - f Q (-k +1)) d^dfg 


By virtue of relation 2.8 the above expression reduces to 


OO 

R y ( Vt 2 ) = J ! V! f + f 0 (l - k)} s x (f, - f 0 k)df, 

K **■ -OO 


With the substitution f^ - kf = f the expression further 
simplifies to 

OO 

R y( fc 1 ^ 2 ) = £ f y*K f + S*f+lf 0 ).t 2 { S x (f)df 

K 1 -00 

3.24 


From eqn. 3.4: 


1 = H°l y ‘ (f + lf o jt 2 } =i 



i27r(f + If )t 

0l H(f + lf 0 ,t 2 )e • 0 
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Because of eqn. 2.17 the left hand side of this equation is the 
response to the input 

i2irft 

A(t)e 3.25 

Denoting this response by Y^(f,t), one obtains finally 
oo 

R y (t v t 2 ) = J Y*(f,t 1 )Y A (f,t 2 )S x (f)df 3.26 

— 00 

This equation has the same form as eqn. 3.8a for stationary 
random input, the only difference being that in case of a modulated 
stationary random input the response Y^(f,t) to A(t)e x is to 
be used instead of the response y(f,t) to e^Trft^ 

The time variable mean square is now 

r° 

<T y 2 ( t ) = Ry(t,t) = J Y*(f,t)Y A (f,t)S x (f)df 3.27 

-00 

and corresponds to eqn. 3.12a. 

In real form equations 3,22 and 3,23 can be used, whereby merely 
y (f,t) and y s (f ,t) must be replaced by the responses Y c ^(f,t); 
Y sA (f,t) to the inputs A(t) cos 2mft and A(t) sin 2irft respectively. 
Prom a computational point of view it is of importance that a 
non- stationary random input of the type considered here can be 
treated in the same way as a stationary random input. 
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4 . Methods of Approximate Solutions 

In Phase I Report the blade angle of attack taken as an 
average over the blade span, had been assumed to represent a 
stationary stochastic process, and the aerodynamic loads on the 
flapping blade had been determined by modulating this stationary 
process with a periodic time function. In a more sophisticated 
theory the blade angle of attack will appear already as a modulation 
of a stationary stochastic process which has to be defined by the 
atmospheric turbulence penetrated by the aircraft at constant 
flight speed. Prom a measurement point of view it is almost a 
necessity to assume the input to be a modulated stationary 
stochastic process, whereby the underlying stationary process 
can be measured with respect to its power-spectral density or 
correlation functions. In contrast, the power- spectral densities 
of non- stationary processes cannot be measured in principle 
and the measurement of their correlation functions requires a 
large set of sample functions, usually not available. The basic 
assumption with respect to the stochastic structure of the input 
made in this report is for lifting rotors from a physical point 
of view plausible, from a measurement point of view almost 
required, and from a mathematical or computational point of view, 
as shown in the preceding sections, a very great simplification. 

The three approximate methods discussed in the following are all 
based on this particular assumption for the stochastic input. 

The presentation Is further limited to a single degree of freedom 
linear system. For the actual lifting rotor more than one 
degree of freedom should be considered, whereby input and 
response would appear in matrix form and where the problem of 
cross correlation functions between the various degrees of 
freedom would occur, as discussed in Phase I Report. Finally, 
the actual lifting rotor description includes non-linearities 
which if small could be considered in the perturbation theory, 
but which would render the general theory presented herein 
inapplicable . 
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4.1 Solution Based on Functional Relation Between Input and 
and Output Double Frequency Spectra 

To illustrate the method, we will consider the approxi- 
mate differential equation for blade flapping, derived in Phase X 
Report ** 

0 + (c 1 + a-|Sint)/§ + (c 2 + a 2 cost)# - (c 3 + a^slntja 4,1 

Here 5 , the blade angle of attack averaged over the blade span, 
is assumed to represent a stationary random process. The right 
hand side of eqn, 4.1 is then a modulated stationary process and 
represents the input to the blade flapping equation. In general 
the factors for 0 (damping of the flapping motion), for 0 (stiffness) 
and for a can also contain other terms of the respective truncated 
Fourier series. 

Taking the Fourier transform on both sides of eqn. 4.1 and denoting 
the Fourier transform of 0 by B, that of by X, one obtains, 
as shown in Phase I Report 

-(2vf 2 ) 2 B(f 2 ) + c^irif^fj) + a { (fg - 5y)B(fg “ 2?> 

" ( f 2 + 2TT^ f 2 + 2?^ f + C 2 B ^ f 2^ + “ 1 B ^ f 2 “ 2r\ + B ^2 + 2 ^ ^ = 

C 3 A ( r 2 ) + — | A ( f 2 “ £7^ + A ^ f 2 + 


Taking the conjugate complex of this equation and substituting 
f 1 for f g, multiplying the two equations and taking the mathematical 
expectation of this product leads to a lengthy functional equation 
between the known power-spectral density for the angle of attacks 


E [A*(f 1 )A(f 2 )J 


0 for f 1 # f 2 
S--(f) for a fg S f 


and the power spectral density for the flapping angle 


E [B*(f 1 )B(f 2 )] =S^ (f 1 ,f 2 ) 


**In later sections constants a t and a Q 

1 it 

respectively and the input modulating f 

a. 

replaced by — g + b sint . 


are replaced by d-j 
unction + a 3 sint 


and 

is 
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For the particular case of eqn. 4.2 the response power spectrum 
S^ (f-pf^) consists only of line masses along the lines f » fg 
and - f 2 * ” ^ ^^Sher order terms of the truncated 

Fourier series in the coefficients of 0 , 0 and a in eqn. 4.2 
are considered, the response spectrum (fpf 2 ) contains also 
line masses along the lines f-j - fg = - , where n = 1,2... 

up to the highest order of the truncated Fourier series. 

Replacing f = kAf, where Af is a small frequency interval 
and k a discrete variable with values k = 1,2,... m, the functional 
equation for S^ (fpf 2 ) can be replaced by a system of complex 
linear equations. Since from physical considerations 

S£ (kAf , lAf ) — * 0 

for sufficiently large k and 1, one obtains a finite number of 
equations to compute the values (kAf, lAf) for all discrete 
values of k,l, for which S^ is assumed to be different from zero. 

The problem is now reduced to the inversion of a large 
number (in the order of many hundred) of linear equations for 
complex variables, having complex coefficients. In Phase I 
Report the relations between the line masses on the diagonal, 

+ 1 

= f 2 and the line masses on the other two lines, f-j - f 2 = - ^ 

had been neglected, and only a relation between the line masses 
on the diagonal line retained. This leads to a real system of 
equations with real unknowns. Since submitting Phase I Report 
it was found that this incomplete system of equations does not 
yield approximation, since not even the trend with increasing 
magnitude of the periodic coefficient ap a 2 , in eqn. 4.1 
is properly established by the incomplete system of equations. 

It was subsequently attempted to solve the complete 
system of equations for the complex unknowns, however for the 
numerical input data considered (see Section 5), the numerical 
experience showed that the complex coefficient matrix is ill 
conditioned and not directly suited for standard iterative tech- 
niques. The program so far completed stores only non zero elements 
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and for any preassigned value of m generates its own complex 
coefficient matrix which is by a standard subroutine split into 
real arithmetic. This representation of a complex matrix in 
real form takes practically twice the original core storage 
but was found to be computationally more convenient and easily 
amenable for double precision. It seems necessary to generate 
a preconditioning matrix, by a trial and error procedure. 

This aspect of the problem, using preconditioning 
matrices, has not yet been exploited. It is unlikely that the 
preconditioning operation and the subsequent inversion of the 
large matrix will be possible without a considerable amount of 
computer time per case. Of physical significance and subject 
to direct measurement is only that portion of the double frequency 
spectrum (fpf 2 ) consisting of a line mass along the diagonal 
f^ = f^, since, according to Section 2.2, only the diagonal line 
mass contributes to the time averaged power-spectral density. 
However, contrary to the assumption made in Phase I, the relations 
involving the line masses on the non-diagonal lines cannot be 
neglected in the computation. 
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4.2 Solution Based on the System Response to Deterministic Inputs 

The response autocorrelation function can 

be either computed from eqn. 3.24, using the response y(f,t) 
to the input e^ 27r ^, or it can be computed from eqn. 3.26, 
using the response Y A (f,t) to the input A(t)e i ^ irft , the actual 
computations to be performed with the equivalent real form 
equations. The latter approach takes less machine time but does 
not economically permit a parametric variation in A(t). For the 
first approach the response calculations are Independent of the 
input modulating function A(t) and variations in this function 
are reflected merely in the double summation of eqn. 3.24. 

In either case a deterministic response analysis over 
a sufficiently wide frequency and time range is required. 

Once the responses are determined, the autocorrelation function 
R y (t^,t 2 ) or the mean square response R y (t,t) requires merely 
a single integration over the applicable frequency range. 

In performing the numerical response computations the 
computer time involved should be an Important factor in selecting 
a suitable one-step or a multi-step method. Truncation and 
numerical Instability problems should not affect the reliability 
of the computations over sufficiently large time intervals. 

It is presumed that the round off errors can be checked with 
a double precision arithmetic - a provision easily available 
in present day computers. 

For second order differential equations as in our blade 
flapping problem. It is possible by a standard substitution 
(Ref. 12, p. 227) to obtain another equation of the same order 
but without the first derivative terms. A multi-step method 
known as Noumerov's method (Ref. 13, p. 137 and Ref. 14, p. 301) 
was used for some sample runs to compute the blade flapping 
response for an input cos27rft, beginning with zero displacement 
and rate of displacement for t = 0. This method has no stability 
problems and the truncation error is of the order of Q(h^), 

Ref. 14, p. 301 ♦ In addition to the known zero displacement 
and displacement rate conditions the method requires one extra 
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start Ing value of the displacement which was computed from 
Taylor series. On an overall basis this multi-step method 
took more machine time than the one -step method described 
below . 

This one-step method was the Runge-Kutta method of 
fourth order specially suited to second order differential 
equations, giving a truncation error of order 0(h) . The 
program is based on the algorithm given in Ref. 12, page 238. 
Being a one-step method it is self starting without any 
stability problems. Numerical comparisons have been made with 
reduced step size and in some other cases with the perturbation 
theory. This rather heuiristic approach toward truncation and 
round off problems indicates that the computational errors are 
too insignificant to affect the reliability of the response 
calculations. 

After numerically computing the deterministic responses 
over an adequate range of frequency and time the quadrature 
operations based on the Simpson rule is carried out in accordance 
with equation 3.24 or 3.26. 

Numerical experience thus far' gained seems to indicate 
that with computer time of 40-45 minutes on machines comparable 
to the IBM 360-50 Model, a reasonably accurate time dependent 
mean square response can be obtained. This assumes that the 
periodically varying damping and stiffness functions are 
explicitly given in the form of Fourier series. Otherwise a 
separate subroutine has to be added to perform such a Fourier 
analysis. 
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4.3 Perturbation Method for Linear Time Variable Systems with 
Small Time Varying Parameters 

A drawback in a complete numerical approach is the 
considerable amount of machine time in computing the responses 
for Inputs cos wt and sinwt over an adequate range of the 
discrete frequency parameter cj . Approximate analytical methods 
on the other hand provide solutions in terms of the variable 
frequency ^ and it is possible to make a qualitative study of 
the response with or without including transient effects. 

To evaluate the correlation function of the response 
by the perturbation method, the solution to known deterministic 
inputs Is expressed as a power series in € , a perturbation 
parameter which in our problem is only a mathematical artifice. 

We have an exact solution when € = 0 but the solution we seek 
is obtained by letting e = 1 . The response of the time variable 
system is then calculated by repeatedly solving the associated 
constant parameter system for known inputs and using the principle 
of superposition. Finally, the computation of the response 
correlation function comprises single quadrature for stationary 
inputs and double quadrature for non stationary Inputs. (See 
equations 3.8 and 2.6) 

A direct Fourier inversion of the correlation function 
to obtain the spectral description of the response, equation 2.4, 
involves computationally inconvenient quadrature operations 
which can be avoided if the time variable parameters and the 
input modulating function are periodic as in the case of a 
flapping blade. In order to make use of the periodicity of the 
system parameters, we express the stochastic response £(t) 
as a power series in e and then obtain the double frequency 
spectrum Sfl according to equation 2.4. This latter 

perturbation scheme is analogous to the one employed for 
non linear stochastic problems. Ref. 2, page 272, except in 
the present blade flapping problem the excitation is a special 
kind of a non stationary process for which the spectral density 
function comprises series of line masses, equation 2.21. 

The response spectrum same as the input spectrum also contains 
line masses. Therefore, for the physically realizable time 
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averaged response spectrum one needs to consider only the 
diagonal terms of the response spectrum. 


4,3,1 Computation of Response Correlation Function 

Consider a time variable parameter linear system 
typified by the equation 


F(0) = 


n n 
0 + Z 

j = 1 


jo, + d, sin( Q q t + 




n H n a 

where 0 = ? ■ ° 

dt n 


4.3 


With the forcing function f(t) for which the spectral density 
function or the correlation function R f (t^ ,t 2 ) is 

known, equation 4.3 takes the form 

F(p) = f(t) 4.4 

Now introduce two more linear operators 

n n n-j 

U0) = 0 + Z c, B 4.5 

j = 1 J 

and 

N(0) = L{0) - F(0) 4.6 

If 

|Cjj » | dj sln( Q o t + g j ) | 4.7 

it is possible to introduce a perturbation parameter e such 
that (10, IT) 

0(t) = 0 o (t) + + e 2 yS 2 (1) + ... 4.8 
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Instead of solving equation 4.3, we now seek a solution to the 
problem 


UP)** f(t) + « N(0) 4,9 

for € = 1 . 

Substituting the power series expansion, equation 4,8, in 
equation 4.9 and equating the coefficients having the same 
powers of € one gets the following infinite system of equations: 


L( 0 O ) = 
L( fiy) = 
L</* k+1 ) 


f(t) 


n 


n-j 


L dj sin( a o t + 9 3 ) 0 O 


n 


n-J 


X dj sin( fl ft t + 9,) 0 


J “ 1 


J 


4.10 

4.11 

4.12 


In equation 4.4, when the random input f(t) is replaced by a 
deterministic harmonic forcing function the response 

y( cj , t) can also be expressed as 


2 

y(w,t) « y Q ( w ,t) +£ y^w^t) + e y 2 («,t) + ... 

The solution of equation 4.10 with f(t) replaced by e iwt gives 


y Q ( (j , t) = H( £j )e 


i<j t 


n 


+ £ B n .(w) e 1 ^ 

J = 1 J 


4.13 


where 



are the roots of the equation 


(IX ) n + £ c ( lX) n-J = 0 
j = 1 3 

and the constants B^(a>) are to be evaluated by satisfying the 
n zero initial conditions of the system. The first part of 
the response In equation 4,13 refers to the steady state 
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solution and the term under the summation sign refers to 
transient solutions which can be neglected for stable systems 
assuming one is only interested in the steady state. If 
transients are of interest, they can always be superimposed 
to the stochastic solution. Henceforth, we will consider only 
steady state solutions. 

When y Q (cj,t) and its n-1 derivatives are substituted 
in equation 4.11 and noting that the right hand side deterministic 
input is periodic, It is possible to solve for y-|(6>,t) in 
closed form. Similarly, y 2 (w,t), y 3 ( W , t) . . .etc . have to be 
solved if correction terms of order more than one are needed. 

As the system Is linear one can set 

y( <j,t) sa. y Q (u ,t) + y«j(cj ,t) +...etc. 4.1 

The correlating function (t^tg) is then obtained from the 
relation 3.8. y («,t) and y (<j,t) in equation 3.20,- re spec - 
tively correspond to inputs coscjt and sintjt or the real and 
imaginary parts of y(w,t) in equation 4.14. 

In the present rigid blade flapping problem the forcing 
function in equation 4.4 Is a separable non stationary process 
of the type of equation 2.16, where the input modulating 
function is periodic. Therefore, the spectral density function 
of the input is given by equation 2.21. In the computation -of 
the response correlation function, the double Integral in 
equation 3.8 reduces to a series of single integrals depending 
upon the number of Fourier terms in the input modulating 
function. Another approach which is better suited for a 
specific problem involving no parametric study of the Input 
modulating function and especially when it is not periodic 
is to compute y A {cj ,t) instead of y(o>,t) and then use relation 
3.26 to compute R^Ct^tgj. The P erturbation scheme to compute 
y A (w,t) is exactly similar to the one described earlier. 

The only difference is that the forcing function now is a 

■j t ^ 

product of the input modulating function and e 
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For any input modulating function A ( t).,y Ao ( w , t) , 
y A1 (<J >t)...y Ak+ i( etc. can be obtained from equations 


Ao 


(w,t) ~f h(r)A(t - r )e iw ^ ” r ^dr 


4.15 


r m 

y A1 (w,t) = / h( r)P 0 (t - r )d t 
— 00 


4.16 


y Ak+1 


OO 

(w,t) =y*h( r )F lc (t -r)dr , etc. 

— OO 


4.17 


where F Q ( t) , . . .F^ t ) represent the right hand side deterministic 
functions in equations 4,11 and 4.12 after replacing the stochastic 
responses /? Q (t) , . . /? k (t) by the deterministic responses y A (<j,t),.. . 

etc# When the input modulating function is periodic, 
convolution integrals in equations 4.15, 4.16 and 4.17 can be 
evaluated in closed form. As the system is linear 

y A ( ^ ,t) ~ y Ao ( w J + y A *| ( ^ ,t) etc . 4.18 

We note in passing that when the input is a stationary 
process, y(w ,t) and y A ( a> ,t) are identical. 

4.3.2 Response Spectral Density Function 

Taking the Fourier transform of equation 4.8 and then 
using relation 2.4, one gets 




( w 


W 2 ) = 


'0 


( w n » o) + € f S 


( + S 


*o*1 U 2 


0 1 0 


( w 


U 


)] 


We are interested only in which for the specific 

type of input considered in this report corresponds to the 
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Physically realizable time averaged power spectrum s l s 
y virtue of relation 2.15, 3.le and 3 fg J ( , \ ^ ' 

expressed as • ( u ) can be 

S /J ( u ) = S 0 ( (o , to ) 

tr “ Sf0 ™ ° f *«“«»» *-10 .t frequency „ glve , 

B 0 U) = H( ^ j 

Application of equation 2.4 gives the spectral den fi1 l - 
in the form P Pral density function 


4.20 


W" 1 ’ w 2 )s f ( 

Similarly the spectral density functions for A(t) .. 
p k +i(t), etc. can be expressed as ' 


4.21 


^k + i (6, i' w 2 ) - Wji , xl^-D^djd^DSn-l-^-K^.^) 


(w 1 - cj n -J( 


O' v “2 " a o ) 


(a) 


n-l 


-1/4 


n n 


J 


W" 1 - fi o’ " 2 -S 0 ) 


J = 1 l£ /- 1 ) n '/d 1 (i)2n-l-J e i( (?J + #i) 
( ^ ^ + 42 ) ( w - ^ n-l 


0 ; v U 2 ~ # 0 ) 


-1/4 


n n 


S ^ k ^1 + a 0 , «2 -fl ) 
J = 1 1 X = ^-^""/^(i^h-l-je-if $ + * ) 


( W 1 - •<? 0 ) n ' J ( " 2 + fl ) n 'l 


S ^ k ( "j ‘ fi o' w 2 +a 0 ) 



+ 1/4 


n 

z 

3 - 1 
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£ (-l) n-J d 1 d 1 (i) 2n-1-J e 1 ^ + 

1=1 J 1 


( u, + + fl.) 


n-1 


S ^ k (w 1 + ^o' W 2 + «o> 

4.22 


Now, take the Fourier transform of equation 4.12 at frequency cj 2 
to yield 


H( U > 2 ) ^k+4 w 2^ 21' ^ W 2 ' ^k^ w 2 " Je J 

-(i) n -J(o> 2 + e 0 ) n - J * k (* 2 + fl 0 ) n - J e i<>J ] 


•jf 

Multiplying both sides with. * k( w -j ) and then taking the 
expectation, equation 2.4, one gets 

S /9k/?k+1 (Ul,W2) = H( " 2 )[ 1 /2 S* k <«v« 2 -« 0 )' j£v e^djU) 1 *-^ 1 

(<*> 2 - fl 0 ) n_J -1/2 S («,, w 2 + i e‘ lej d j (l) n ' J ' 1 (a> 2 +fl o ) n ' J J 

4.23 


With the help of relations 4.21, 4.22 and 4.23, the response 
spectral density function can thus be obtained to any desired 
order. 
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5. Numerical Examples 

The equation studied "by the perturbation and numerical 
methods reads 

0 + (c 1 + d^ sin Q 0 t) j? + (c 2 + d 2 cos Qq t) 0 ~ A(t)a(t) 5.1 

Selecting a time unit for which the rotor angular velocity^ = 1 
and substituting c-j = ^ , , c 2 = 2c -j , d 2 = d-j and 

A(t) = c-j + 2d 1 sint, an approximate blade flapping equation 
valid up to moderate advance ratios is obtained. (For details 
of equation 5.1 refer to Phase I Report, page 17). For practical 
rotors the non dimensional-inertia number y varies from 2-10, 
therefore we have assumed a typical value of y — 4 in the numerical 
examples. Note also that the system parameters c 1 and d 1 
which are linearly related to the advance ration also appear 
in A(t) . However, the system parameter d^ is varied from 0 to 1 
only in the left hand side of equation 5.1 without changing 
the right hand side deterministic function A(t). The computational 
scheme with different values of d^ is thus associated with two 
specific functions 


and 


A(t) = 0.5 +0.4 sint 


2 

a 0.25 + W 


As the function A(t) which corresponds to the actual physical 
system at an advance ratio of 0.3, is not simultaneously changed 
along with the system parameter d 1 , the computed values of the 
time variable mean square response (t,t) and the time averaged 
response spectra (<j) correspond to the actual blade flapping 
problem only for d 1 = 0.2. The intent of this report is not 
so much to carry out an extensive parametric study of y and M 
for different input spectra which are of interest in the atmospheric 
turbulence study but to establish the range of validity of the 
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perturbation method with regard to fi and then within this 
admissible range develop a computational scheme to compute the 
time variable mean square response and the time averaged spectra 
of the flapping oscillations. With the computer program written 


in Fortran IV language, it is possible to evaluate R^ (t,t) 
for any desired values of y and V when the stochastic input 


a(t) has a known spectral density function Sg (<j ) . 

Now, coming to the actual description of the computational 


scheme and presentation of computer results, it is convenient 


to describe them in three stages: 


In the first stage, the range of validity of the perturbation 
method with respect to the advance ratio fi has been established 
by comparing the perturbation system responses with that of the 
Runga-Kutta method results. The computer program, for any 
preassigned values of c-j and d 1 , gives the system responses 
according to these two methods provided the input to the 

cl 14- 

system is of the form (-§ + "b-j sint)e iw ^. Numerical results 

presented in Figures la to 1c refer to three typical frequency 

values of 0.5, 1 and 1 .5 with constants a =2 and b, =0. 

o l 

A comparison of system responses. Figures la to lc, Indicates 
that up to d^ « 0.7 (or approximately - 1) the perturbation 
scheme should provide reasonably accurate results. 

In the second stage, the time variable mean square 
response R^(t,t) is evaluated according to equation 3.27, 
by integrating ’the product of the spectral density function 
&=* ((*>) and the square of the absolute value of the system 
response to the input (— g + b^ sint)e . The limits of 
integration in equation 3.2 have been truncated to -3 to 


+3. This finite frequency range of integration seems to be 
adequate for applied purposes because in the numerical examples 
discussed here the value of the integrand for |oj| > 3 is less 

_3 , 

than 10 . . The quadrature routine is based on Simpson* s rule 

with a stepsize of 0.1. Figure 2 refers to two cases - 


stationary and non stationary inputs to a constant parameter 
system. The first case corresponds to the blade flapping 
equation at zero' advance ratio and as expected, the mean 



square response, dotted lines in Figure 2, are time invariant. 

The second case, though not directly relevant to the present 
blade flapping problem is an analytical model of considerable 
interest where the mean square response is time dependent due to 
non stationarity in the input. The strong time dependency of 
the mean square response, full lines in Figure 2, is due to the 
fact ghat the time invariant or the constant part of the input 
with — ^ = 0.5, is not small compared to 0.4 sint, the time 
variable part. For a spectral description of such a separable 
non stationary process refer to Section 2.3, Figure 3 shows the 
time dependency of the mean square response of time variable 
parameter systems subject to non stationary excitations. Here 
both the time variability of the system parameters and non 
stationarity of the input contribute toward the non stationarity 
of the response. 

The third stage comprises the computation of the time 
averaged response spectra according to equations 4.21, 4.22, 

4.23 and 4.20, Figure 4 summarizes these numerical results 
for different values of the system parameter d-j . For two 
extreme values of d^ =0 and d^ = 0.8, Figure 5 shows the 
comparison between the perturbation values and the NASA 
conducted simulator results, and Figure 6 also refers to a 
similar comparison with d-j = 0.2. Considering the discrepancy 
in the time averaged input power spectrum between the simulator 
study and exact analytical values. Figure 5, and also other 
types of errors due to finite filter band width, etc. inherent 
to simulator results, the perturbation values within the admissible 
range of the perturbation scheme mentioned earlier agree reasonably 
well with the simulator results. Figures 5 and 6. 



Response 


0 + (0.5 + d sint) £ + (1. + d cost) fi ~ cost 




FIGURES .la-1 ct 



(For negligible difference between the Runga-Kutta and the perturbation method 
values the responses are shown by full lines only. Otherwise, the dotted lines 
represent the Runga-Kutta method values and the full lines perturbation method 
values.) 
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Figure 2: Mean Square Response of a Constant Parameter System Subject to Stationary and 

Non Stationary Excitations 
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Figure 3: Study of the Effects of Time Varying Parameters on the Mean Square Response. 

(The case d = 0.2 refers to the blade flapping problem at an advance ratio fi - 0.3 
and inertia number y = 4) 
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6 • Conclusions 

Same as in Phase I Report it is assumed here that 
atmospheric turbulence produces an input stochastic loading of 
the lifting rotor blades in forward flight which is of the 
separable kind and consists of a stationary random process 
modulated by a periodic time function. Since submitting 
Phase I Report the problem of random blade flapping has been 
further studied and the following new results have been 
obtained: 

6.1 The approximate method of solving the functional relation 
between input and output double frequency power spectral 
densities for the flapping rotor blade in forward flight, 
tentatively suggested in Phase I Report, has been checked 
against simulator results. While the method gives the 
correct order of magnitude effects of moderate advance 
ratio on the response power spectral density, it cannot 
be used as a quantitative estimate of these effects. 

6.2 In the approximate method of 6.1 all off-diagonal terms 
of the double frequency input and output power spectral 
densities were neglected. A more elaborate approximation 
was tried, including the off-diagonal line masses. The 
resulting system of linear equations has a complex coeffi- 
cient matrix which turned out to be ill conditioned and 
not directly suited for iterative techniques. 

6.3 Solutions have been ' deve loped based on the system response 
to deterministic inputs. Such responses can be computed 
with standard numerical methods like the Runge-Kutta 
method. The deterministic input consists of a harmonic 
forcing function with or without being modulated by the 
right - hand side deterministic time function. If the time 
variability of the system parameters compared to the asso- 
ciated constant parameters is not too large, the deterministic 
response can be obtained with adequate accuracy also by 

the perturbation method. In either case the response 
autocorrelation function or the time variable mean square 



^43- 


response can be computed by a single integration over the 
applicable frequency range. 

6.4 A perturbation method has been developed in two forms. 
According to the first form the deterministic responses 
are computed by adding to the solution for the associated 
constant parameter system corrections from the time varying 
parameters. The deterministic responses thus obtained 

are then inserted into the appropriate integrals over the 
applicable frequency range, representing autocorrelation 
function or time variable mean square. 

According to the second form of the perturbation method 
the response power spectral density for the associated 
constant parameter system is first computed and then improved 
by adding the necessary corrective cross spectral terms. 

While the first form of the perturbation method is best 
suited for the computation of the time variable mean square 
response, the second form of the perturbation method lends 
itself particularly well to the computation of the time 
averaged response power spectral density. 

6.5 The numerical examples presented herein have the main 
purpose to determine the range of applicability of the 
perturbation method- and to evaluate for some typical 
assumed cases the stochastic structure of the blade 
flapping response. The results of this method are compared 
to the results of NASA conducted simulator studies. Also 
some typical response time histories obtained with the 
perturbation method are compared to those obtained with 
the more elaborate Runge-Kutta numerical integration 
method. On the basis of the numerical examples treated 

it can be concluded that the perturbation method of 
determining blade responses to stochastic inputs is 
approximately valid up to an advance ratio of one in 
combination with a Lock inertia number of 4. 
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So far the studies were primarily concerned with 
questions of methodology in treating time variable 
systems under non-stat ionary stochastic inputs of the 
separable type. The methods have been applied to a 
simplified approximate differential equation of blade 
flapping in forward flight. No effort was made to 
predict the stochastic input from a given atmospheric 
turbulence structure. Probably • such a prediction will 
require empirical parameters to be obtained from model 
or flight tests. The extension of the studies to multi- 
degree of freedom representations of the blades will 
require data on cross correlation functions between the 
generalized stochastic loads, which also should be based 
on experiments. A rather straightforward extension of 
the present studies concerns a more accurate representation 
of the blade in flapping or flap-bending with more complex 
expressions for the time variable damping and stiffness 
terms. Such an extension together with the computation 
of typical random response data over a wide range of blade 
parameters and stochastic inputs is presently in work. 


\ 
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